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Abstract. 

The hemispheric pattern of solar filaments is considered in the context of the global 
magnetic field of the solar corona. In recent work Mackay and van Ballegooijen have shown 
how, for a pair of interacting magnetic bipoles, the observed chirality pattern could be 
explained by the dominant range of bipole tilt angles and helicity in each hemisphere. 
This study aims to test this earlier result through a direct comparison between theory 
and observations, using newly-developed simulations of the actual surface and 3D coronal 
magnetic fields over a 6-month period, on a global scale. 

In this paper we consider two key components of the study; firstly the observations 
of filament chirality for the sample of 255 filaments, and secondly our new simulations of 
the large-scale surface magnetic field. Based on a flux-transport model, these will be used 
as the lower boundary condition for the future 3D coronal simulations. Our technique 
difi'ers signiflcantly from those of other authors, where the coronal field is either assumed 
to be purely potential, or has to be reset back to potential every 27 days in order that 
the photospheric field remain accurate. In our case we ensure accuracy by the insertion of 
newly- emerging bipolar active regions, based on observed photospheric synoptic magne- 
tograms. The large-scale surface field is shown to remain accurate over the 6-month period, 
without any resetting. This new technique will enable future simulations to consider the 
long-term build-up and transport of helicity and shear in the coronal magnetic field, over 
many months or years. 
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1. Introduction 

Solar filaments are large regions of dense cool plasma suspended in the hot 
corona (Engvold, 1998). They are well observed in Ha, either in absorption 
against the disk or in emission at the limb. For historical reasons in the 
latter case they are usually termed "prominences". Coronal magnetic fields 
are key to the existence of filaments (Mackay, 2005), as they both support the 
filament mass against gravity and help to insulate it from the surrounding 
corona. 

Since the earliest magnetograph observations (Babcock and Babcock, 
1955) it has been known that filaments form above polarity inversion lines 
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Figure 1, Definition of filament chirality, showing the axial magnetic field (B) direction, 
and barb orientation, in a dextral filament and a sinistral filament. The magnetic polarities 
in the photosphere beneath are also indicated. 



(PILs) in the photosphere, which divide regions of positive and negative 
magnetic flux. In the filament itself, the dominant magnetic field component 
is found to be along the filament's long axis (Foukal, 1971; Leroy, 1989). 
The filament may be classified as having either dextral or sinistral chirality 
(handedness) according to the direction of this dominant axial field (Mar- 
tin, Bilimoria, and Tracadas, 1994). Viewing the filament from the positive 
polarity side of the PIL, a right-pointing field is dextral and a left-pointing 
field is sinistral. This is illustrated in Figure 1. 

The property of filament chirality is notable because it follows a hemi- 
spheric pattern, whereby dextral filaments dominate in the Northern hemi- 
sphere and sinistral in the Southern hemisphere. The pattern was initially 
observed in polar crown prominences by Rust (1967), who noted that the 
field direction was opposite to that expected from differential rotation. It was 
later observed by Leroy, Bommier, and Sahal-Brechot (1983) with Hanle- 
effect measurements, and then confirmed by Martin, Bilimoria, and Tracadas 
(1994) and Pevtsov, Balasubramaniam, and Rogers (2003) using Ha obser- 
vations. The latter paper found that 80% - 85% of quiescent filaments (which 
form between the remnants of decaying active regions) followed the rule. A 
feature of the hemispheric pattern is that it remains the same when the 
Sun's polar field reverses, approximately every 11 years. 

Zirker et al. (1997) suggested that the hemispheric pattern could be 
produced by a combination of supergranular motions, differential rotation, 
and magnetic reconnection (Martens and Zwaan, 2001), without recourse to 
unobservable subsurface phenomena as in some earlier theories (van Balle- 
gooijen and Martens, 1990; Rust and Kumar, 1994; Priest, van Ballegooijen, 
and Mackay, 1996). This suggests that the pattern for quiescent filaments 
could be a natural result of the observed evolution of bipolar active regions, 
and their interaction as they are transported poleward. A series of papers, 
beginning with van Ballegooijen, Cartlcdgc, and Priest (1998), have used 
coupled numerical modelling of the surface and coronal magnetic fields to 
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investigate this idea. Having reached the conclusion that emerging active- 
region hehcity does play an important role (Mackay and Gaizauskas, 2003), 
Mackay and van Ballegooijen (2005) gave a convincing explanation for the 
hemispheric filament pattern. However, this was a parameter study which 
considered only a pair of model magnetic bipoles in a localised region of 
the solar corona. In this study we aim to extend the model to the global 
corona, and to simulate the real Sun for a period of 6 months, allowing us to 
compare the local direction of the simulated magnetic field with the chirality 
of a large number of observed filaments during the period. 

A distinction must be made between two approaches that have been used 
to model the coronal magnetic field based on photospheric flux transport. In 
the first approach the photospheric flux transport simulations are carried out 
independently from the coronal field modelling {e.g. Mackay and Lockwood, 
2002; Wang, Shcclcy, and Lean, 2002; Schrijvcr and Derosa, 2003). Here the 
flux transport simulations are used at discrete intervals, perhaps every 27 
days, to supply the lower boundary condition for the extrapolation of a po- 
tential magnetic field in the coronal volume. The potential field is analysed at 
that time instant and then discarded, a new potential field being constructed 
after the next 27 days. In this approach, which has been successful in the 
past, there is evidently no direct coupling between the flux transport model 
evolving the photospheric field, and the extrapolated coronal fields. This 
precludes any memory in the coronal field of previous field line connectivity. 

In contrast, more recent flux transport models have coupled the photo- 
spheric and coronal magnetic fields, and include the models of van Ballegooi- 
jen, Cartledge, and Priest (1998) and van Ballegooijen, Priest, and Mackay 
(2000), developed to investigate the hemispheric pattern of filaments. The 
technique described in the latter paper has been used with small modifica- 
tions in a series of papers (Mackay, Gaizauskas, and van Ballegooijen, 2000; 
Mackay and van Ballegooijen, 2001, 2005, 2006), and is used again in this 
study. As carried out in Mackay, Gaizauskas, and van Ballegooijen (2000), 
the radial magnetic field component in the photosphere is taken as the 
initial boundary condition, and an initial potential coronal field extrapolated 
from it. Both the photospheric and coronal fields are then evolved together, 
allowing the build-up of magnetic energy and helicity in the coronal volume 
result of the photospheric motions. 

Mackay, Gaizauskas, and van Ballegooijen (2000) applied the model to 10 
filaments forming above an observed activity complex. In these simulations, 
the emergence of new flux was not included, so the photospheric field became 
less accurate over time. For this reason the photospheric fields were reset to 
the observed field after every 27 days of evolution. At this point, the sheared 
coronal field had to be discarded and a new potential field constructed. 
Although a significant advance on previous models with a purely potential 
coronal field, the build-up of magnetic energy and shear was limited to rather 
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short periods of time compared to that on the real Sun. Despite the loss of 
accuracy in the surface field, the authors tried running the simulation for 
several months without resetting, and found in fact that the chirality results 
were slightly improved. 

We have now modified the model, enabling us to simulate coronal mag- 
netic fields for much longer periods of time without having to reset the 
photospheric or coronal fields every 27 days. Key to our new technique is 
the insertion of individual active regions throughout the simulation, based 
on the observed photospheric field. The new surface flux transport model 
is described and illustrated in this paper, simulating the real Sun over a 
particular 6-month period in 1999. 

A sample of 255 filaments have been observed in Ha images from Big Bear 
Solar Observatory (BBSO), over the same period. The chiralities of these 
filaments have been determined where possible, as described in this paper. 
In a follow-up paper the full 3D coronal field evolution will be computed, 
coupled to the new surface simulation. The resulting magnetic skew will be 
compared with the chirality of filaments at the locations where they were 
observed. 

The structure of this paper is as follows. Section 2 describes the filament 
data, and the observed chiralities. The numerical model for surface flux 
transport is outlined in Section 3, while Section 4 describes in detail the 
method for including newly-emerging active regions. Example surface simu- 
lation results are given in Section 5, and a concluding summary in Section 
6. 



2. Filament Observations 

The filament data set consists of 255 filaments identified in Ha images 
from Big Bear Solar Observatory, over a period from May to October 1999. 
These are large, stable filaments, traditionally classified as quiescent or in- 
termediate (Engvold, 1998), rather than unstable short-lived "active region 
filaments". However, while the majority of the filaments are located between 
active regions, our selection criteria do include a small number of filaments 
forming within active regions. All filaments in the sample were observed for 
at least 4 consecutive days, but the majority had lifetimes much longer than 
that. 

The filament locations have been overlayed on synoptic magnetograms 
from the US National Solar Observatory, Kitt Peak, for the corresponding 
Carrington rotations CR1949 to CR1954. For example. Figure 2 shows the 
filament locations for CR1952 and CR1954. Thick white lines indicate the 
filament locations, and the background shading shows radial magnetic field 
Br on the solar surface, with horizontal coordinate longitude and vertical 
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Figure 2. Observed filament locations overlayed on Kitt Peak synoptic magnetograms for 
two selected Carrington rotations, CR1952 and CR1954. On the magnetograms white 
represents positive flux and black negative flux. 



coordinate sin (latitude). The filament locations span a range of latitudes, 
from the equatorial regions to the polar crowns (60° - 70° latitude). Notice 
how all of the filaments overlie PILs in the photospheric magnetic field. The 
chiralitics of the filaments have been determined by the technique described 
in Section 2.1, using the Ha images. 
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Figure 3. Example sequence of Ha observations for a single filament, taken on four con- 
secutive days from BBSO full-disk Hq images. White squares highlight barb locations 
identified on the images. From the barb directions this filament is clearly seen to be 
sinistral. 

2.1. Filament Chirality 

Each individual filament was followed for up to 7 daily observations about 
central-meridian passage. An example sequence is shown in Figure 3 for one 
particular filament. On each day the angle of view is altered by the Sun's 
rotation, and a number of barbs may be observed. The orientation of these 
barbs with respect to the main axis of the filament may be used to determine 
the filament's chirality, since filaments with right-bearing barbs are dextral, 
and those with left-bearing barbs are sinistral (Martin, 1998). 

The filament in Figure 3 is clearly sinistral as all of the observed barbs, 
shown by white squares, are left-bearing. However, the classification is not 
always so straightforward, as some filaments appear to have barbs of both 
chiralitics (Pevtsov, Balasubramaniam, and Rogers, 2003). This may be a 
result of short-lived perturbations in the local filament structure, or may 
simply be due to the viewing angle as the filament rotates with the Sun. 

To minimise subjectivity in the classification, the following statistical 
technique was used. For each filament, the number of dextral barbs and 
the number ns of sinistral barbs was determined, over all available Ha images 
for that filament. A t-test (details in Appendix A) was used to determine 
whether the filament chirality was significant, based on the observed number 
of barbs of each type. 
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Figure 4- The 255 observed filaments and their chirality on a time-latitude plot. Squares 
are dextral filaments and asterisks are sinistral filaments. Small dots indicate filaments 
where the chirality could not be unambiguously determined from the available Ha data. 



With a critical value oi t = 1.5 in the statistical test, 99 of the 255 
observed filaments could be classified as either dextral or sinistal from the 
available observations. Of these 64 were dextral and 35 sinistral. Investi- 
gation revealed a further 24 filaments to have identifiable chirality, but 
not enough barbs to pass the t-test at this level. Thus the full data set 
contains 123 filaments of known chirality and 132 of unknown chirality. The 
classification of each observed filament is shown in Figure 4, along with its 
latitude and date of observation. 

Notice that there are a number of exceptions to the hemispheric pattern. 
Among the filaments of known chirality, 88.7% follow the hemispheric pat- 
tern in the Northern hemisphere, while in the Southern hemisphere 73.1% 
follow the pattern. This result is in line with those of Martin, Bilimoria, 
and Tracadas (1994) and Pevtsov, Balasubramaniam, and Rogers (2003). 
The aim of this project is to compare directly the chirality of filaments 
in our sample with the simulated coronal magnetic field. To allow such a 
comparison we need to model the photospheric field accurately over the 
observed period, and this is described next. 
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3. Surface Model 

Accurate simulation of the photospheric magnetic field is a key part to the 
project, providing the lower boundary condition that drives the 3D coronal 
field in our model. In this section the flux transport model is described, 
while the technique for including newly-emerging flux is given in Section 4. 

3.1. Flux Transport Equations 

The transport of magnetic flux on the surface of the Sun has been well 
studied over a number of years and a standard model has emerged (Sheeley, 
2005 and references therein) . Concentrated flux regions originating at active 
latitudes are spread over the solar surface by convective motions and sheared 
by differential rotation. It is now well accepted that there is also a poleward 
bulk flow, or meridional circulation. These three effects may be incorporated 
in an advcction-diffusion equation for the normal magnetic field on the solar 
surface (Leighton, 1964; Devore, Boris, and Sheeley, 1984; Wang, Nash, and 
Sheeley, 1989). In this paper we write the magnetic field as B = V x A and 
solve for the vector potential A. The evolution of Bj- on the solar surface is 
then specified in spherical polar coordinates (r, 6,(f)) by the following pair of 
equations: 

dAe _ D dBr 

dt ~ '''^ ' Re sine d<t> ' ^ ' 

dt ^ Rq ae ^ ^ 

Here D is a diffusion constant that models the random walk of magnetic flux 

due to the changing supergranular convection pattern (Leighton, 1964; an 
alternative is the flux-dispersal model of Schrijver, 2001). For our simulations 
a value of D = 450km^s~^ was used. The differential rotation velocity = 
fl{0)RQsm0 is well determined from observations. We use the Snodgrass 
(1983) formula for the angular velocity which, in the Carrington frame, 
gives 

J^(^) = 0.18 -2.3 cos^^- 1.62 cos'' ^ deg day"^ (3) 

The meridional flow is an order of magnitude weaker than the differential 

rotation so it is observationally less well known (Hathaway, 2005). As in 
Mackay and van Ballegooijen (2006), we use the proflle 



Ug{9) = C COS 



71" (^max + ^min ^ 
2 (6'max - ^'min) 



cose, (4) 



where C is chosen to give a maximum flow at mid-latitudes of 16 ms ^. 
In addition to the advection and diffusion terms, a source term is added 
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to Equations (1) and (2) to represent newly-emerging active regions. This 
allows us to simulate the photospheric evolution continuously for many 
months. 

3.2. Numerical Method 

To compute the evolution of Bj. on the solar surface the flux transport 
equations for Aq and A^f, are solved numerically in a two dimensional domain 
covering both hemispheres, where 10° < 9 < 170° and 0° < (j) < 360°. For 
ease of computation, transformed coordinates {x,y) are used, where 

x = (f>, y = - In (tan (0/2)). (5) 

The number of grid points in each direction is (361, 281), giving a resolution 
of 1° at the equator. At each stage terms on the right-hand sides of Equations 
(1) and (2) are computed using finitc-diffcrcnccs, before advancing in time 
with a time step dt = 60s. Spatial derivatives are obtained to second-order 
accuracy by defining variables on a staggered grid, where A and j are defined 
on cell ribs and B on the cell faces. 

Boundary conditions are applied to Br and are simply periodic in the 
x-direction and closed in the y-direction. The initial condition Br{x,y) is 
derived from observed magnetograms as described in Section 3.3. 

3.3. Magnetogram Data and Initial Surface Distribution 

As observational input to the simulation we use readily available synoptic 
magnetogram data from the US National Solar Observatory, Kitt Peak. 
The resolution of this data is sufficient for the present study as we are 
interested in the large-scale global magnetic field, and in quiescent filaments 
which may span tens of degrees on the Sun. For the surface simulation the 
magnetograms are used in two ways: firstly to obtain an initial condition for 
the simulation, and secondly to determine where new large-scale flux has 
emerged (Section 4). 

Each normal component magnetogram is taken in daily strips over a 27 
day period (a single Carrington rotation). This creates a problem when a 
simultaneous map is needed at all longitudes, as is required for the initial 
condition in the global simulation. McCloughan and Durrant (2002) deal 
with the problem by using a synoptic evolution equation, rather than the 
flux transport equation itself. We take an alternative approach and produce 
an initial configuration directly from the magnetogram. A simultaneous 
magnetic map of the photosphere is created by correcting the observed 
magnetogram for differential rotation, as illustrated in Figure 5 which shows 
the creation of the initial condition for our simulation. The effects of diffusion 
and meridional flow are neglected as they are much weaker over a timescale of 
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(a) yiei 




Otl947 CR1948 GR1949 




Figure 5. Correction of synoptic magnetograms for differential rotation, producing the 
initial condition for the simulation, starting on day 106. Three consecutive maps (a) are 
used to produce a single corrected map (b), applying the formula (6). The dashed line 
4> = 0j.g£ shows the reference longitude for day 106, and remains unrotated. The dot-dashed 
lines show the effect of the rotation on the boundaries of the original CR1948 map. 



13 days, the time between either the left or right edges of the magnetogram 
and the central reference longitude ^^-gf. This is the longitude viewed at 
central meridian on the start day of our simulation. 

The required correction is dependent on longitude in the magnetogram, 
which effectively corresponds to time over the 27-day period. In Figure 5 we 
have chosen the reference longitude (/>j.g£ = 180°, which was observed mid- 
way through the 27 days (corresponding to day of year 106 in this case). 
If the whole map is to represent this single day, then points to the right of 
(/)j.gf must have extra rotation applied, as their central-meridian passage was 
earlier than 0j.gf . Accordingly, points to the left of ^j,g£ must have rotation 
removed. If the original point is at longitude and colatitude 9, then the 
new longitude in our simultaneous map will be 

(pniO) = ^ref + (^o - ^^-ref)- (6) 

Here Q{6) is the differential rotation profile (3) in the Carrington frame, 
and r^o = 13.2 deg day~^ is the rotation rate of the Carrington frame with 
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respect to an observer on Earth. Notice that all points remain at the same 
colatitude during this operation. The amount of rotation depends on their 
colatitude, as is clearly seen by the change in shape of the original map 
boundaries in Figure 5. It can also be clearly seen in Figure 5 that we must 
also make use of the two neighbouring magnetograms to produce a single 
map using this operation. Thus for our initial condition based on the mid- 
way point of Carrington rotation CR1948, we must use magnetograms for 
CR1947, CR1948, and CR1949. 

While this technique only gives an approximate map of what would have 
been on the surface on the start day of our simulation, the simulation results 
(Section 5) show that it is sufficiently accurate to reproduce the surface field 
with future evolution. 



4. Emerging Flux 

In this section we detail the method for selecting newly-cmcrgcd active 
regions and inserting them into the simulation. Over the 6-month timescale 
of our simulation, a significant amount of magnetic fiux will be added by new 
active regions, driven by the solar dynamo (Charbonneau, 2005). Schrijver 
and Harvey (1994) estimate a rate of flux emergence from bipolar regions 
of 6.2 X lO^^Mxday"^ near solar maximum. They find that enough flux 
emerges in these active regions to replenish the entire photospheric flux in 
about 4 months. For our 6-month simulation, we must clearly include newly 
emerging active regions in order to reproduce the surface magnetic field, and 
subsequently the coronal field, with any realism. 

Solar active regions are thought to result from the buoyant rise of sub- 
photospheric magnetic flux tubes which break through the surface in an ft 
loop, producing the characteristic bipolar magnetic configuration (Parker, 
1955; Babcock, 1961). These magnetic bipoles are oriented in an approx- 
imately cast-west direction, although the bipolc axis is typically tilted so 
that the leading polarity is nearer to the equator (Hale et al, 1919; Wang 
and Sheeley, 1989). The vast majority of bipoles in the Northern hemisphere 
have the same leading polarity, and those in the Southern hemisphere have 
the opposite polarity (Hale and Nicholson, 1925). With the global reversal 
of the Sun's polar magnetic field every 11 years, all the polarities reverse. 

Since the model does not include a description of the magnetic field inside 
the Sun, we cannot predict when and where magnetic flux will emerge. As 
the ultimate aim is to reproduce observed coronal features, it would be 
inappropriate to insert new active regions at random locations and times 
during the evolution, as has been used in some previous surface flux trans- 
port models (van Ballegooijen, Cartlcdge, and Priest, 1998; Mackay and 
Lockwood, 2002; Baumann et al, 2004; Mackay et al., 2004). Instead the 
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observed magnctograms throughout the simulation period are analysed to 
determine where and when new flux has emerged (Sheeley, Devore, and 
Boris, 1985). Magnetic flux is then inserted into the simulation in the form 
of magnetic bipoles with properties matching those of the observed active 
regions. Of course, such a technique based on the synoptic magnetograms 
will only pick up the larger emerging flux regions, but we believe that this is 
sufficient to produce the large-scale field over long time-scales, as is required 
in this study. 

4.1. New Flux Regions 

Locations of new flux emergence are found by comparing successive synoptic 
magnetograms. Some care is necessary in doing so because the field will have 
evolved in the intervening 27 days. The comparison is made easier by apply- 
ing the correct amount of differential rotation to the earlier magnetogram, 
as in Section 3.3. Figure 6 illustrates the procedure for a particular case. The 
earlier observed magnetogram (CR1948) is shown in Figure 6(a), and it has 
been corrected for differential rotation with reference longitude ^j,g£ = 0° 
(corresponding to the final day of that rotation). The resulting magnetogram 
is shown in Figure 6(b). This is the magnetogram which is compared with 
the observed magnetogram for the next rotation (CR1949), shown in Figure 
6(c). In Figures 6(a) to (c), white areas represent magnetic flux above 50 G 
and black areas represent flux below — 50G. 

A semi-automated procedure has been developed to determine the loca- 
tions of new flux, based on some simple criteria. This is found effective in 
reducing the time taken for a manual search. The stages in this procedure 
are as follows: 

1. Determine the boundary coordinates of all regions of flux on the later 
observed magnetogram. These regions are deflned as continuous pixel 
groups with radial fleld \Bj.\ greater than a threshold of 50G. 

2. Compute an absolute difference map between the later observed magne- 
togram and the rotated earlier magnetogram. An example is shown in 
Figure 6(d). Here red areas represent flux (of either sign) present only 
in the later observed map, and blue areas represent flux present only in 
the rotated earlier map. 

3. Select all regions from flrst stage where there is a majority of red flux 
on the difference map. These represent new regions. 

4. Flag any anomalous regions for manual attention. These include regions 
where there is mostly only one polarity of flux, regions that do not appear 
to follow the Hale polarity law, and also any regions not considered to 
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Figure 6. Stages in processing the synoptic magnetograms for the detection of new bipolar 
regions: (a) observed map for CR1948, (b) result of applying 27.27 days of differential 
rotation, (c) observed map for CR1949, and (d) absolute difference map. fn (a), (b) and 

(c) , white represents magnetic flux above 50 G and black represents flux below — 50G. fn 

(d) red indicates new areas of fiux and blue indicates old areas. 

be new which have very high peak strength. All probable new regions 
which have not been flagged are then classified as "definitely new". 

5. The flagged regions are included or excluded manually, and the user can 
select any further regions as desired, in order to improve the accuracy 
of the simulation. 

The final step has been particularly useful in the case of large activity 
complexes, which may consist of multiple bipolar regions in close proximity. 
On some occasions it has been necessary to add in a new bipole in existing 
active regions where more fiux appears to have emerged. A total of 119 
bipolar regions were determined in the 6 Carrington rotations CR1949 to 
CR1954 inclusive. Of these, 97 were detected by the automated program, 
the remainder being added manually. 

Having located the bipolar regions where new fiux has emerged, the 
following properties of these regions are measured: day of central-meridian 
passage, coordinates of bipole centre, tilt angle, half-separation distance be- 
tween peaks, and flux. The locations of the bipole peaks were determined by 
taking the centroid of fiux of each polarity. Figures 7(a) and (b) summarize 
the data for the 119 observed new regions. Figure 7(a) shows the typical 
relation between bipole tilt angle and latitude of emergence known as Joy's 
Law (Hale et al., 1919). We find that the average tilt angle, from a least- 
squares flt, varies as 0.38 x latitude, but with a large spread of data about 
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Figure 7. Summaxy of properties of the 119 bipolar regions: (a) tilt angle 6 against bipole 
centre latitude, and (b) magnetic flux against half-separation distance p between centres 
of each polarity. In both plots the dashed line is a linear least-squares fit. 



the regression line as seen in Figure 7(a). This is in good agreement with 
Wang and Sheeley (1989), who found a slope of about 0.4, and a similarly 
large spread. In Figure 7(b) we see the expected positive correlation between 
total magnetic flux and active region size (Schrijver and Harvey, 1994). 

Our 119 emerging regions contribute a total flux |$_|_| + |$_| of 1.18 x 
10^^ Mx over the 177-day simulation. This agrees with the estimate of Schri- 
jver and Harvey (1994) for flux emergence from bipolar regions near solar 
maximum. It is significantly higher than the total flux present in the photo- 
sphere at any one time, which at the start of our simulation is 5.09 x 10^^ Mx. 
Emerging flux is clearly vital to the photospheric field evolution over this 
6-month period, and must clearly be included in our simulations. 

4.2. Bipole Description 

Emerging flux is included in our simulation by the insertion of magnetic 
bipoles with properties matching those of the observed regions. The mag- 
netic bipoles all take the same functional form, as used in Mackay and 
van Ballegooijen (2001) and Mackay and van Ballcgooijen (2006). This is 
a three-dimensional field, allowing not only the insertion of the bipole at 
the photospheric level, but also in the corona. In the corona this bipole field 
may be non-potential, representing the helicity contained in newly-emerging 
regions on the Sun. In our computational (x, y, z) coordinates a bipole is 
given by the vector potential 

A, = pBoe'^-'ze-^^, (7) 
Ay = Soe°-Ve-^, (8) 



arysurface.tex; 1/02/2008; 16:09; p. 14 



Modelling the Global Solar Corona 15 

A, = -pBoe^-^xe-''^, (9) 

where 

^=i-' + ^y^ + y\ (10) 

and the parameter p is the half-separation distance between the peaks of the 
two polarities. This is for a bipole with zero tilt angle; for non-zero tilt angle 
6 the field is simply rotated. The twist parameter /3 controls the helicity 
of the bipole. In the present paper only the photospheric footprint of this 
field will be required, which does not depend on /?, but the helicity will be 
important in the future coronal simulations. 

For the present we set z = and disregard the Az component. This gives 
the radial surface magnetic field 



P 



+ 

p2 



(11) 



which is illustrated in Figure 8(a) for a non-zero tilt angle 6 corresponding to 
a rotation of the coordinates (.x, y) in Equation (11). Note that in Equation 
(11) Br is independent of the twist parameter /?, so that changing (3 alters the 
coronal field but gives the same photospheric field. Finally, the observations 
give the total flux $ of the bipole rather than the peak field Bq, but for the 
model bipole the two are related by 



^° = 7^7^- ^^^^ 



4.3. Insertion of New Regions 

New bipoles are inserted by adding their corresponding vector potential, as 
described in the previous section. There are two complicating factors; firstly 

the exact day of emergence is unknown, as some bipoles may emerge on the 
unobserved side of the Sun, and secondly there may be pre-existing flux at 
the emergence location. 

Because the observed synoptic magnetograms have a time cadence of only 
27 days, there is considerable uncertainty about the exact emergence date 
of new flux regions. In particular it is unlikely that a region appeared on 
the exact day it was observed at central-meridian, and it may well have 
emerged on the unobserved far-side of the Sun. In view of this, and to allow 
time for coronal structures to develop before the date of filament observation 
at central-meridian passage, we have emerged all bipoles 7 days before their 
date of observation. To ensure accuracy we have "evolved" the parameters 
of each bipole back in time by 7 days (see Appendix B for more details). 
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Figure 8. Illustration (a) of the bipole given by Equation (11), with tilt angle 5 and 
half-separation p. An example insertion from the simulation is seen "before and after" on 
(b) day 250 and (c) day 251. All three plots show magnetic field strength Br, with solid 
contours being positive flux and dashed contours negative flux. 

Seven days is an arbitrary choice, but we have varied the number of days 
and found no critical effect. This aspect of our model could be improved in 
future, either by more accurate modelling of the development of individual 
active regions, or by using newly-developed helioseismic far-side imaging 
techniques (Lindsey and Braun, 2000). Note however that such techniques 
are not yet sufficiently developed for use in this present work. 

Now consider how a new bipole will be inserted in a region with a 
considerable pre-existing magnetic field. High coronal field lines previously 
connected to the photosphere in that region find themselves with nowhere 
to connect down after the new bipole is inserted. To avoid any associated 
problems in the three-dimensional simulation we have adopted a technique of 
"sweeping" to move any earlier strong field out of the region of insertion. This 
is not without some physical basis, as one might expect a newly emerging 
flux tube to distort the previous coronal field out of its path (Yokoyama and 
Shibata, 1996; Krall et al, 1998). 
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While the sweeping takes place the simulation time is frozen, so that 
the bipole insertion is still instantaneous. The ordinary flows are turned 
off, and a new outward velocity of the form Vg = exp{—As) is applied in 
the photospheric and coronal region surrounding the bipole centre, where 
the coordinate s represents distance from the bipole centre. The constant 
A is chosen so as to give zero flow at the edge of the insertion region. This 
sweeping velocity advects the field out of the insertion region, and sweeping 
stops once the field near the bipole centre reaches a suitably low level. We 
have found that a threshold of O.OSBq works well, where Bq is the maximum 
magnetic field of the new bipole. 

Figures 8(b) and (c) show the insertion of a new region on day 251 of the 
simulation. There is a considerable pre-existing field on day 250 (Figure 8b) , 
and the bipole emerges in what was largely a region of positive flux (solid 
contours), but with an encroaching region of negative flux (dashed contours) 
at higher latitude. Considerable sweeping was needed to reduce the field 
strength in the insertion region. This has resulted in the pre-existing negative 
flux being pushed up to higher latitude to make way for the new bipole 
in Figure 8(c). The bipole insertion procedure in 3D will be described and 
illustrated when the 3D coronal field simulations are discussed in a following 
paper. This will demonstrate how our newly inserted bipoles reconnect with 
the overlying coronal field to produce a plausible magnetic configuration. 

By the method outlined in this section we have successfully incorporated 
the insertion of newly-emerging bipoles into the photospheric simulation, 
allowing continuous evolution over many months. Furthermore, the bipoles 
are three-dimensional and non-potential in nature. This represents a signif- 
icant advance over previous models, and in the subsequent 3D simulations 
will allow us to study the build-up of magnetic helicity and shear in the 
coronal field over long periods, driven by the emergence and evolution of 
active regions. 



5. Surface Simulation Results 

In this section we illustrate the newly developed surface simulations over a 
particular 6 month period from April to October 1999. This will form the 
continuous photospheric boundary condition for the global coronal mag- 
netic field simulations to be described in a following paper. The simulation 
starts from day-of-year 106, mid-way through Carrington rotation CR1948. 
It then runs continuously for 177 days until day 283 at the end of CR1954. 
Two different simulation runs are illustrated: one where no emerging flux 
is inserted during the simulation, and one where the 119 bipolar regions 
described in Section 4.1 are inserted at the relevant days. In both cases the 
initial condition is the corrected magnetogram for day 106 (Figure 5b). 
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Figures 9 to 11 illustrate the results at different stages in the simula- 
tion. In each figure, (a) is the observed synoptic magnetogram, (b) is the 
simulation run with no bipolc insertions, and (c) is the simulation run with 
bipole insertions. All three of these figures correspond to the final day of 
a particular Carrington rotation, rotations CR1949, CR1951, and CR1954. 
The observed magnetogram has been corrected for differential rotation (Sec- 
tion 3.3) with reference longitude ^^g£ = 0°, giving an approximation to 
the instantaneous field of the Sun on the final day of the given Carrington 
rotation. This is to facilitate comparison with the simulated field, and the 
observed magnetogram has also been smoothed for easier comparison. In 
each plot the horizontal coordinate is longitude and the vertical coordinate 
is sin(latitude). All plots use the same levels for Br, with a saturation level 
of ±50 G. Black areas represents negative flux and white areas positive flux. 
For each day, the zero contour of Br from the observed magnetogram (a) is 
overlayed on the simulated magnetogram (c) for comparison. 

After 1 month of evolution the simulation with no bipole insertions (Fig- 
ure 9b) reproduces the diffuse background field well, particularly nearer 
to the poles. This magnetic fiux results from the dispersal of earlier active 
regions which have been carried poleward by the meridional circulation. The 
relatively slow speed of the circulation, at most 16ms~^, ensures that flux 
emerged in the active region belts since the start of the simulation has not 
yet reached these higher latitudes. At lower latitudes there are numerous 
regions of strong fiux missing from this simulation; these are the active 
regions that have emerged since the start of the run. 

When emerging regions are included the simulation reproduces the ob- 
served field at day 147 well, as evidenced by visual inspection and the 
overlayed zero contour on Figure 9(c). The extra bipolar regions in the right- 
most 90° of the simulated magnetogram (to the right of the black dashed 
line) should be ignored in the comparison. These belong to the observed 
magnetogram for the next rotation, CR1950, and so do not appear in Figure 
9(a) which is the observed magnetogram for CR1949. They appear on this 
day in the simulation because all of our new regions are inserted 7 days 
before their date of observation at central-meridian passage. 

By day 202 (Figure 10) the simulation with no bipole insertions still repro- 
duces the field well near the poles, but it has lost all stronger flux regions at 
lower latitudes. However, the simulation with emerging flux maintains good 
accuracy. It is less accurate in complex flux regions, such as in the top right 
of Figure 10(c), but the large-scale picture is still reasonable. In principle the 
structure in these regions could be improved by carefully inserting a larger 
number of small bipoles, but this would be labour intensive. In any case, 
such regions become smoothed out over time and the small-scale complexity 
is lost. 
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Figure 9. Comparison of simulation results at day 147: (a) observed magnetogram cor- 
rected for differential rotation, (b) simulation with no bipole insertions, and (c) simulation 
with bipole insertions. In all three plots, white indicates positive flux and black negative, 
with a saturation level of ±50 G. On the simulated magnetogram (c), the zero contour 
from the observed magnetogram (a) is overlayed in white. To the right of the black dashed 
line, new regions from the next rotation will appear in (c) but not in (a), because we insert 
them 7 days before observation at central meridian. 
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Figure 10. Comparison of simulation results at day 202: (a) observed magnetogram cor- 
rected for differential rotation, (b) simulation with no bipole insertions, and (c) simulation 
with bipole insertions. In all three plots, white indicates positive flux and black negative, 
with a saturation level of ±50 G. On the simulated magnetogram (c), the zero contour 
from the observed magnetogram (a) is overlayed in white. To the right of the black dashed 
line, new regions from the next rotation will appear in (c) but not in (a), because we insert 
them 7 days before observation at central meridian. 
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Figure 11. Comparison of simulation results at the end of the simulation, on day 283: (a) 
observed magnetogram corrected for differential rotation, (b) simulation with no bipole 
insertions, and (c) simulation with bipole insertions. In all three plots, white indicates 
positive flux and black negative, with a saturation level of ±50 G. On the simulated 
magnetogram (c), the zero contour from the observed magnetogram (a) is overlayed in 
white. 
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The contrast between the two simulations is even more evident on day 
283, where the simulation with no bipole insertions (Figure lib) is now 
beginning to lose accuracy even near to the poles, with the field very weak. 
When bipoles are inserted (Figure 11c) the simulation still gives a good 
reproduction of the large-scale field, containing both strong newly emerged 
flux and the diffuse remnants of earlier active regions. Therefore it can be 
seen that our technique for reproducing the surface field is able to maintain 
accuracy to the observed field. 

6. Conclusion 

In this paper we have developed a new technique for simulating the observed 
photospheric magnetic field on a global scale over many months, without 
having to reset to observed magnetograms every 27 days. In a follow-up 
paper these new simulations will be coupled with simulations of the 3D 
coronal magnetic field, where our motivation is to explain the hemispheric 
pattern of the axial magnetic-field direction in solar filaments, a property 
known as filament chirality. 

The essential features of the new photospheric magnetic-field simulations 
are as follows: 

1. A surface flux transport model, based on synoptic magnetogram data 
from the real Sun. 

2. Accuracy is maintained not by resetting to the observed field, but by the 
emergence of new active regions throughout the course of the simulation. 

3. These active regions are in the form of magnetic bipoles, with proper- 
ties chosen to match those of bipolar regions observed in the synoptic 
magnetograms. 

We do not reset the photospheric field because in the future 3D sim- 
ulations this would remove memory of previous field-line connectivity in 
the coronal field. It would prevent the long-term build-up of non-potential 
magnetic energy and helicity, which are key to the formation of filament 
magnetic structures. An advantage of the magnetic bipoles we insert is their 
mathematical form, which allows them to be inserted also into the corona in 
the 3D simulation. They can then be given a non-zero magnetic helicity, and 
such active-region helicity is thought to be important for filament chirality 
(Mackay and Gaizauskas, 2003; Mackay and van Ballegooijen, 2005). 

Our semi-automated technique for finding new bipolar regions has been 
found to work well, with the program detecting automatically 97 of the 
119 regions inserted in the 6-month simulation. The extra regions have been 
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chosen manually; this has largely been necessary for areas with many bipoles 
in close proximity, such as in activity complexes. With the 119 emerging 
regions, the simulation is successful in reproducing the large-scale photo- 
spheric field, even over a period of 6 months. The accuracy of the simulation 
is seen clearly when the observed polarity inversion line is overlayed on 
the simulated magnetograni. Larger areas of positive and negative magnetic 
flux are well reproduced, as is the shape of the inversion line at most of the 
observed filament locations. 

It is these locations where our simulated 3D magnetic field will be com- 
pared with the chirality of observed filaments. In a recent paper, Mackay and 
van Ballegooijen (2005) have shown how the observed properties of active 
regions on the Sun, along with the observed large-scale surface motions, 
can produce the observed hemispheric pattern. However, they considered a 
simple pair of bipolar magnetic regions, simulating only a localised region of 
the solar photospheric and coronal magnetic field. In the present study, we 
will test their conclusion by comparing global simulations with real filament 
data, over a 6-month period. In the global simulations, the 3D coronal field 
will be coupled directly to the continuous evolution of the photospheric field, 
so the project relies heavily on the accuracy with which we have been able 
to simulate the observed photospheric field in this paper. 

For comparison with the simulation, we will use a filament data set which 
includes up to 7 daily observations of each of 255 filaments, taken from BBSO 
Ha images over the 6-month period. As described in this paper, we have 
determined the chiralities of individual observed filaments using a statistical 
technique based on the orientations of their barbs. The chirality has been 
unambiguously determined as either dextral or sinistral for 123 filaments, 
and it is at these locations where we will compare the observed chirality 
with the direction of skew in the 3D coronal magnetic field simulations. 
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Appendix 



A. Statistical Test for Filament Chirality 

A statistical t-test is used to determine filament chirality from classification 
of individual barbs. Suppose a single filament has n classifiable barbs, giving 
n observations Xi,X2, ■ ■ ■ ,Xn, where Xi = +1 or —1 according as the barb 
is right-bearing or left-bearing respectively. Then the number of dextral 
barbs is = Xj, and the number of sinistral barbs is = n — 

n^. It is reasonable to assume that follows a binomial distribution with 
parameters {n,p), and we assume p = 0.5. 
The test is based on the statistic 

f- ^d-"P (Al) 
^/np{l - p) ' 

which should be near to if neither chirality is significant. The classification 
scheme is then 

t > T dextral, 
t < — T sinistral, 
|t| < T not classified, 

where T is some chosen threshold. For large enough n, t should approximate 
a normal distribution with mean n and variance 1. As our sample sizes are 
relatively small we choose T = 1.5. For a filament with 9 barbs, Equation 
(Al) shows that this corresponds to a threshold ratio n^/n = 0.75. 



B. Backward Evolution of Bipole Parameters 

To insert a bipole 7 days before its date of observation at central-meridian 
passage, the bipole properties must be "evolved" back in time for 7 days, 
so as to match the observations as closely as possible when it does reach 
central meridian. The properties in question are the location of the bipole 
centre (^cen, (pcen), the tilt angle 5, the half-separation p between centres of 
polarity, and the magnetic fiux. 

The movement of the centres of each polarity may be inferred from 
the formula (3) for differential rotation 0(0), as meridional flow may be 
effectively neglected over such a timescale. This allows the evolution of the 
bipole centre, tilt angle, and half-separation to be calculated. Suppose the 
bipole properties are known at time t = to- Then if the centre of the leading 
polarity is at {6i,<pi{t)), and that of the following polarity is at {02,<p2{t)), 
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it can be shown that 

<Pcen(t) = (pcenCio) H 2 

^cen(t) = ^cen(to) for all t, (B2) 

2p(to)cos5(to) + i?0smecenf^o*' 
= \l f + ' i^l sin^ gcen + (p(to) sin ^(to))^ 



i?0 sin 6'cen 2 



(B4) 



where 



Qo = ^^(^+)-J^(^-), (B5) 
and the colatitudes 0+ and 0- are given by 

e^ = e,,^ + PM^^ and ^_ = ^cen - ^^^^^1^. (B6) 

The magnetic flux above the 50 G threshold will also change, and depends 
on both advection and diffusion. It has been determined approximately by 
test simulations running in the forward-time direction. By varying the lat- 
itude, tilt-angle, and half-separation, a look-up table was created relating 
the initial and final flux for different parameters. It was not necessary to 
include different values of initial flux in this look-up table, as varying the 
initial flux simply results in a proportional scaling of the flux at all times. 
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